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In the present work, we propose a new set of coherent structures that arise in nonlinear dynamical 
lattices with more than one components, namely interlaced solitons. These are waveforms in which 
in the relevant anti-continuum limit, i.e. when the sites are uncoupled, one component has support 
where the other component does not. We illustrate systematically how one can combine dynamically 
stable unary patterns to create ones such for the binary case of two-components. In the one- 
dimensional setting, we provide also a detailed theoretical analysis of the existence and stability 
of these waveforms, while in higher dimensions, where such analytical computations are far more 
involved, we resort to corresponding numerical computations. Lastly, we perform direct numerical 
O , simulations to showcase how these structures break up, when exponentially or oscillatorily unstable, 

2^ ' to structures with a smaller number of participating sites. 
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C/^ ! INTRODUCTION 

' One of the highly active areas of investigation of Hamiltonian nonlinear systems over the past decade has been 
the examination of nonlinear dynamical lattices of the discrete nonlinear Schrodinger (DNLS) type. Chiefly, this 
development has arisen due to the multitude of applications of pertinent models that have emerged in areas such as 
nonlinear optics and atomic physics. 

More specifically, in the optical context, the setting of fabricated AlGaAs waveguide arrays [H has been one of 
the most prototypical ones for the application of DNLS models. There, the interplay of discreteness and nonlinearity 
revealed many interesting features including Peierls-Nabarro potential barriers, diffraction and diffraction management 
0, and gap solitons 0], among others; see also the reviews [3, [1] and references therein. 

Another recent development, which also promoted the analysis of discrete systems in connection with nonlinear 
optics was the proposal p and creation 0, |8|] of optically induced photonic lattices in photorefractive crystals such 
as SBN. This paved the way for the observation of a large set of exciting nonlinear wave related phenomena in such 
crystals. As a representative subset, we mention here the formation of patterns such as dipole quadrupole (TT 
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■ and necklace llj solitary waves, impurity modes [l^ , discrete vortices 13|, , rotary waves [l5| , higher order Bloch 

^ ' modes [3| and gap vortices [l3|, two-dimensional (2D) Bloch oscillations and Landau-Zener tunneling 18|, wave 

•'-j formation in honeycomb ^1^, hexagonal [20| and quasi-crystalline lattices and recently the study of Anderson 

rS localization in disordered photonic lattices |22| . Although this setting is mostly studied in the continuum context with 

. a periodic potential (and sometimes in the presence of the inherent crystal anisotrow), it has also spurred a number 

■ ■ ■ ' of studies in the DNLS context with the saturable photorefractive nonlinearity 23l. |24|. 



Lastly, another physical realization of such nonlinear dynamical lattices arose over the past few years in atomic 
physics through the examination of Bose-Einstein condensates (BECs) trapped in periodic potentials. There, once 
again, a reduction of the relevant model can be formulated in the tight-binding approximation within the mean-field 
limit, reducing the so-called Gross-Pitaevskii equation with a periodic potential to a genuinely discrete nonlinear 
Schrodinger equation [2^. 

In both the nonlinear optical and in the atomic physics setting discussed above, multi-component systems were 
also examined in recent investigations. More specifically, the first observations of discrete vector solitons in optical 
waveguide arrays were reported in [2^, the emergence of multipole patterns in vector photorefractive crystals was 
presented in [2711. w hile numerous experiments with BECs were directed towards studies of mixtures of different spin 
states of ^'^Rb~[28l. [2§| or ^^Na [30| and even ones of different atomic species such as '^^K-^'^Rb [31j] and "^Li-^^^Cs 
(3^ . It should be noted that while the above BEC experiments did not include the presence of an optical lattice, the 
addition of such an external optical potential is certainly feasible within the present experimental capabilities [33[ . 

Our aim in the present work is to propose and analyze a family of solutions particular to multicomponent (in 
particular, binary, although more-component generalizations are certainly possible) systems of DNLS equations. We 
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dub these proposed solutions "interlaced" discrete solitons and vortices, a name stemming from the feature that the 
profiles of the modes in the two interacting components will have a vanishing intersection of excited sites in the 
extreme discrete limit of zero coupling between adjacent nodes of the lattice. In these structures, the first component 
will be excited where the second component is not and vice-versa. In the one-dimensional case, we show how to 
interlace in a stable fashion simple, as well as more elaborate, bound states of the system 3J]. For such solutions, we 
consider their existence and stability properties also from an analytical point of view, using as a starting point the 
anti-continuum limit (of no-coupling between the sites). Then we generalize our considerations to higher dimensional 
settings, showcasing the potentially stable interlacing of more elaborate structures, such as discrete vortices [^] (but 
also of vortices with non- vortical structures). We present detailed stability diagrams of such interlaced structures, 
and also examine their dynamics when they are found to be unstable. 

Our presentation is structured as follows. In section II, we present the model and general mathematical setup. In 
section III, we illustrate both analytically and numerically the properties of such structures in Id settings. In section 

IV, we generalize these considerations to a numerical investigation of higher dimensional settings. Finally, in section 

V, we summarize our findings and present our conclusions. 



MODEL EQUATIONS AND MATHEMATICAL SETUP 

We consider a set of coupled DNLS equations 



0, 
0, 



(1) 



where n is a Z?-Dimensional index and is the discrete Laplacian in D dimensions. We look for stationary solutions 
{it„}, {vn} through the relation 



Un{t) = exp(iAit)u„, Vn{t) = exp(jA2i)?;„. 
The dynamical equations ([1]) then transform into 

- Aiu„ + (gii|w„p +gi2|w„nu„ + CA£)M„ = 0, 

-A2W„ + (.gi2|w„|^ +.g22|Wn|^)Wn +C'A_DW„ = 0. 



(2) 



(3) 



The stability is determined in a frame rotating with frequency Ai for Un{t) and A2 for Vn(t), i.e., we suppose that 
Ur, (t) = exp(iAit) [w„ + ei') (t)] , K (t) = exp(i A2i) [i;„ + (t)] . (4) 
The small perturbations ci'°''(i), with k = 1, 2, can be expressed as 



^n \t) — a„exp(iAi) -I- 6„ exp(— jA*t), ^n Hi) — c„exp(iAt) -f d„ exp(— jA*t), 
leading to the linear stability equations 

with 
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FIG. 1: (a) Profiles of |01) interlaced solitons with gi2 = 0.5 and C = 0.15. (b) Dependence on C of the real and imaginary 
parts of eigenfrequencies of small perturbations about |01 > with gi2 — 0.5. Dashed lines correspond to Lyapunov-Schmidt 
predictions of equations (|20p . and (|2ip . (c) Two-parameter stability diagram in the plane of intersite (C) and inter-component 
((?12) coupling, indicating regions of occurance of Hopf bifurcations (H), exponential instability (E), and stability domain (S). 



Soliton and vortex solutions are calculated using methods based on the anti-continuous limit. Upon calculating 
these solutions at C = 0, we continue them to finite coupling by varying C or other parameters (such as the interspecies 
nonlinearity strength 512). 

We are interested in interlaced solitons (ISs) in ID lattices and interlaced vortices (IVs) in 2D and 3D lattices. The 
excited sites at C = are equal to u and v, except for a phase factor exp(i^), while UnVn = at the corresponding 
excited site. These values are 

M = 0, VAi/ffii, C = 0, ^^2/522- (9) 

In what follows, we choose Ai = A2 = A and gu = 522 1- We also choose gi2 < 1 as, for 512 > 1 interlaced 
solitons and vortices are unstable for every value of C. 

ANALYTICAL AND NUMERICAL RESULTS FOR ID INTERLACED SOLITONS 

Existence and stability 

We consider interlaced solitons which are labeled by \AB >= \A > \B >, where A,B = 0, 1, 2, . . .. This number 
indicates the "order" of the excited state at the anti-continuous limit, whose phase ^ = 0,7r is chosen so that the 
isolated solitons (i.e. when gi2 = 0) are stable for any small C. For instance, the ground state |0 > means m„ = u6n,o 
and the first excited state |1 > will be taken to mean u„ = u{Sn,i — Sn.-i) at the AC limit. Thus, the state |01 > 
corresponds to u„ = udn,o, = v{Sn,i - Sn-i) and |12 > to ii„ = u{Sn^i - Sn.-i), v„ = v{dn,2 + - 6n,o- 

We first analyze the |01 > state, which is stable for C < Cq. At C = Cq the ISs become unstable through Hopf 
bifurcations (the value of Co differs as a function of the rest of the system parameters such as gi2, however the above 
scenario is robust). Cascades of this type of bifurcations arise as C increases and, when, C > Ci, the ISs become also 
exponentially unstable. There is a special region for 912 G [0.27, 0.37] where the system experiences an inverse Hopf 
bifurcation recovering the stability in a window. The system becomes unstable again through Hopf bifurcations for 
gi2 S [0.27,0.34] and exponential instabilities for gi2 S [0.35,0.37]. Besides, for 512 G [0.38,0.47] there exist windows 
with only exponential instabilities. Fig. [1] illustrates all of the above features, by showcasing a typical example of the 
[01 > state, a typical continuation of its principal linearization eigenfrequencies A, and a full two-parameter diagram 
of the stability of this state in the two-parameter plane (C, 312). 

For [12 > states, the scenario is essentially similar to the [01 > case, although, in essence, it is considerably simpler 
due to the absence of any inverse Hopf bifurcations and restabilization windows. Fig. [2] shows the corresponding 
features for [12 >, as Fig. [T]for the [01 > case. 

Dynamics of unstable solitons 

First, we analyze the dynamics of [01 > ISs. Fig. [3] shows the evolution of a typically unstable (i.e. oscillatory 
unstable) [01 > IS with g = 0.2 and C = 0.6. The oscillatory evolution of the instability eventually transforms the 
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FIG. 2: (a) Profiles and (b) dependence on C of the real and imaginary parts of eigenfrequencies of small perturbations of |12) 
showing the same features and for the same parameters as in Fig[T] Dashed lines correspond to Lyapunov-Schmidt predictions 
of equations (|24|l . and (|2ip . (c) Two-parameter stability diagram in the plane of intersite (C) and inter-component (312) 
coupling. 




FIG. 3: Time evolution of the density of the two components for a slightly perturbed unstable [01 > IS with gi2 = 0.2 and 
C = 0.6. 



mode into a |00 > state, which is a stable state of the system. The final excited site is typically the same for the 
{Un} and {Vn} coordinates, although in some cases (even for the same parameters set), the asymptotic excited site 
does not need to be same. However, the amplitude of for the nth site is not identical, i.e. 7^ In a similar 
vein, Fig. 2] shows the evolution of an oscillatory unstable |12 > IS with g — 0.2 and C = 0.4, and, analogously to the 
|01 > case, the IS evolves to a |00 > state (although the finally populated site is not the central one of the original 
configuration). 




FIG. 4: Time evolution of the density of the two components for a slightly perturbed unstable |12 > IS with gi2 = 0.2 and 
C = 0.4. 
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Perturbation analysis 



In this subsection, we attempt to understand in some more details the above observed results of the numerical 
computations in connection to the stability properties of the interlaced soliton solutions. More specifically, we evaluate 
explicit expressions of the interlaced solitons' eigenvalues for the configurations discussed above. The method is based 
on the expansion in the coupling constant C , in the vicinity of the anti-continuum limit. 

In the limit C = 0, as illustrated above, there are two types of solutions, i.e. Un — Vn — 0, and the non-zero solutions 
given by Eqs. ([9]). In this limit, one can also easily notice that the eigenvalue problem ^ will give 

A = ±A, ±A(1 - .gi2/.gii), ±A(1 - 512/522) (10) 

for the zero solutions and 



A = ±0 



(11) 



for the non-zero solutions Q. 

It can be directly inferred from the analysis of the underlying linear problem that the stable eigenvalues A = ±A will 
expand creating a band of continuous spectrum when C is increased. Therefore, this eigenvalue will not be discussed 
further. The instability for a soliton solution will then be determined by the bifurcation of the remaining eigenvalues. 

Let us now first consider the profile of |01 > ISs. It is clear that for finite C the solutions will be deformed from 
their AC-limit profile. The leading-order solution up to 0(C) is then found to be 



"0 a/ — 



Ui = u_i = 



^Agii(l-gi2/322) ' 



Vo =0, 1)1= = 



\/Ag22 ' 



(12) 



The next step is to consider the stability problem when the coupling is turned on. To the leading order, the 
eigenvalue problem of this particular configuration is then given by 



(13) 



where 



a = diag(J), 
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(14) 



and IdixA is the identity matrix of size 4x4. 

Since we have expanded u„ and Vn in a power series of C, then it is natural that we also expand all the involved quan- 
tities in C, i.e. M = Mo + CMi+C'^M2 + 0{C^), E = Eo + CEi+C'^E2+0{C^) and A = Aq + CAi + C^Aa + O^C^). 
It can be checked that Mo is a singular self-adjoint matrix. 

Substituting the expansions to the eigenvalue problems (jl3p will give us to the leading order 



Mo Eo = Ao cr^o, 



(15) 



from which one will obtain that Aq is given by Eqs. (fTO|) and (fTTj) . In the following, let us first consider the case of 
Ao = which are of three pairs, with the corresponding eigenvalues of Mq Eq = denoted by ej, j — 1, 2, 3. Therefore, 
one can write 



So = ^ Ci e 



The next order equation of (|13p gives us 



MoEi = Aicr^io --A^iSq. 



(16) 
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Using the Fredholm alternative theorem, the above equation will have a solution if the right hand side is orthogonal 
to the null space of A^Oi which it is. Hence, the value of the correction Ai cannot be obtained yet and a solution Si 
of Hn]) can therefore be calculated for any Ai. 

The equation of order 0{C^) from ^ can be easily deduced to be 



Mo ^2 = A2 cr ^0 + Ai cr ^1 - Ml ^1 - M2 ^0 



(17) 



Projecting the equation above to Cj, j = 1, 2, 3, i.e. basis of the null space of Mq, will give us the following eigenvalue 
matrix 



-2gi 







-2gi 



\ (3ii-9i2)A 

which can be immediately solve to yield 



(3ll-9l2)A (Sll-S12)A 



-2gii Q -2gii 

(Sll-Sl2)A 



Ai = ±0, ±0, ±2 



Cl 

C2 I = - 
C3 




ffll - 312 

This illustrates that there is a pair of eigenvalues bifurcating from zero as given by 



A = ±2C 



511 



511 - .912 



(18) 



(19) 



(20) 



The same procedure can be applied to bifurcations of the non-zero eigenvalues. In this case, the calculation is even 
simpler as applying the Fredholm alternative to the 0{C) equation of (|13p already gives us a solvability condition 
from which we obtain that bifurcating eigenvalues are 



A = ±(1 - gi2/ff22)(A + 2C), ±(1 - gi2/.gii)(A + 2C), 



(21) 



The above procedure can also be similarly and immediately applied to the configuration |12 > ISs. The only 
difference is that for that solution one will obtain a stability matrix M of size 28 x 28. 
For 1 12 >ISs, we can obtain the solution in a power series of C as 



liO =0, Ul = — 

Wo = -V2 = -V-2 = - 
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(22) 



Continuing to finding the eigenvalues, we will also immediately obtain that in place of (|18p . one will obtain the 
following eigenvalue problem 
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from which we can obtain eigenvalues bifurcating from zero as 



(3ll-3l2)A 



-2311 
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(23) 



A = ±. 



1 - 512/51 
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1 - 512/511 



1 -512/5: 



-C. 



22 



(24) 



Bifurcations from the non-zero eigenvalues for this case can also be shown to yield Eq. (PT|) . 

The above analytical expressions give us a detailed handle on the dependence of the relevant eigenalues on the 
system parameters. Comparisons of the analytical results obtained here with the numerical ones are presented in 
Figs. [T][2] where one can see that the analytical expressions are in relatively good agreement with the numerical 
results. It should be noted that although such analytical considerations are procedurally straightforward to generalize 
in higher dimensions, the relevant calculations are extremely tedious and will thus not be pursued here. Instead, 
we now turn to numerical computations to showcase the existence and potential stability of interlaced solitons and 
vortices in higher-dimensional settings. 
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(a) (b) (c) 




FIG. 5: (a) Plot of real (top) and imaginary (bottom) parts of interlaced vortices with gi2 ~ 0.5 and C — 0.1. (b) Dependence 
on C of the real and imaginary parts of eigenfrequencies of small perturbations about such solutions with gi2 — 0.5 (b) 
Two-parameter stability diagram in the plane of intersite (C) and inter-component (</i2) coupling. 



(a) (b) (c) 




FIG. 6: (a) Plot of real (top) and imaginary (bottom) parts and (b) dependence on C of the real and imaginary parts of 
eigenfrequencies of small perturbations about interlaced soliton-vortex solutions, showing the same features and for the same 
parameters as in Fig. [5] (c) The corresponding two-parameter stability diagram. 



NUMERICAL RESULTS FOR INTERLACED STRUCTURES IN HIGHER DIMENSIONS 

For the case of 2D lattices, we consider two different interlaced structures. On the one hand, we examine interlaced 
vortices (IVs) whose configurations in the AC limit are given by mo,i = u, uo,-i = ui o = iu, M-i,o = —iu; 
vi^i — ~iv, = "D, = iv, ~ —V. On the other hand, we also study a discrete soliton interlaced with 

a vortex (IVSs) whose configurations in the AC limit are given by mo,i = u, uo,-i = ~u, uifi = iu, w_i,o — —iu, 

IVs experience a set of bifurcation scenaria which are qualitatively similar to those of the 1 12 > ISs. IVSs experience 
the same scenario as well, with the basic difi'erence that they appear to exist for all C's (within the range examined 
i.e., up to C = 2) for g < 0.4. Also, notably, the IVSs experience solely Hopf bifurcations in a fairly small region 
inside the exponential-|-Hopf region. Figs. [5] and [H] summarize the corresponding findings in a way similar as for the 
Id configurations, presenting not only typical profiles of the modes, but also typical mono-parametric continuations, 
as well as their full two-parameter stability diagram in the space of inter-site and inter-species coupling. 

In the case of 3D lattices, we consider two interlaced vortices conjoined in the shape of a cube. In the AC limit, this 
cube is given by = u, = —u, — iu,, = ~iu; — v, — —v, = iv, 

V-i^-i^i — ~iv. The structure is stable near the AC limit, with the size of the window of stability diminishing as gi2 
approaches 1, and with instability setting in via Hopf bifurcations. In the 2-parameter continuation figure shown in 
Fig. El the couphng is only continued to C = 0.75, but it is observed that for values between gi2 ~ 0.703 and gi2 = 1, 
the instability further degenerates into Hopf and exponential instabilities. It should also be noted that within the 
region of instability, there some exist isolated points or very narrow regions where inverse Hopf bifurcations may be 
observed, which have been omitted from the graph for clarity. Let us note in passing here that the interlaced vortices 
in the "vortex cube" shown in Fig. [7] are perhaps not the prototypical interlaced structure that one would expect in 
3D; instead one might expect a structure where each vortex is confined in a diagonal plane within the cube (with the 
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FIG. 7: (a) The top set of four panels show an interlaced cube {g — 0.85) in a grid of size 11 x 11 x 11 that has been continued 
to coupling C = 0.20. The level contours shown correspond to Re(u„,m_i) = Ke{vn,m,i) ~ ±0.5 max in blue and 
red (dark gray and gray, in the black-and-white version) respectively, while the imaginary ones, lm{u„^rn,i) = Im(w„^m_j) = 
±0.5 max {tt„,m,i}, are shown by green and yellow (light and very light gray, in the black-and-white version) respectively. The 
bottom panel shows the real eigenfrequencies of small perturbation, (b) The top panel shows the stability diagram, while the 
bottom panel shows the imaginary eigenfrequencies of small perturbation. 



two such planes intersecting transversally). We were, however, unable to trace such a structure even in the vicinity 
of the anti-continuum limit. 



Dynamics of unstable structures 

The dynamics of the oscillatory unstable IVs in 2D lattices with gi2 = 0.2 and C = 0.3 is shown if Fig. [51 The 
evolution results in the transformation of the original structure into single-peaked or multi-peaked solitons. Excited 
peaks do not coincide for J7„ and Vn- The vorticity of each vortex is lost. Fig. [9]shows the dynamics of an oscillatorily 
unstable interlaced vortex-soliton structure with gi2 — 0.2 and C = 0.45. This mode evolves spontaneously towards 
single-peaked solitons. The excited peaks are in the same site in both lattices in this example. 

Dynamics of the interlaced cube with 512 = 0.85 is shown in Fig. [TOl Here the coupling is continued to C = 0.6. 
This is well past the threshold of stability for this value of gi2, and takes the configuration into the region of both 
exponential and oscillatory instabilites. It is observed that when a peturbation of magnitude 0.01 is applied, only a 
single site survives for long times (in this case for Vn). 

Although these are prototypical results of the dynamical evolution, which we have generically observed to lead to 
less elaborate (and often purely single-peaked) structures in this setting, it should be stressed that the specific details 
of the unstable dynamical evolution of each structure depend considerably on the values of the parameters, as well as 
partially on the type/strength of the perturbation. 



CONCLUSIONS AND FUTURE CHALLENGES 



In the present work, we have illustrated the possibility to successfully interlace structures which are stable in each one 
of the components (either simple ones, such as single site solitary waves, or more elaborate ones, such as bound states 
and vortices) in order to produce stable multi-component interlaced solitons/vortices. We have continued the resulting 
structures from the anti-continuum limit of no inter-site coupling to finite coupling and illustrated the intervals of 
stability, as well as the ones of both exponential and oscillatory (Hopf) instabilities. We have given detailed two- 
parameter diagrams of the stable ranges of the solutions as a function of the inter-site and inter-component couplings. 
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FIG. 9: Snapshots showing (a) \U„{t)f and (b) for unstable interlaced vortex-solitons with gi2 — 0.2 and C — 0.45. 
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These revealed that the linear stability of the interlaced structures necessitates sufficiently weak coupling (typically 
no larger than 0.4, with the relevant range decreasing as the inter-component interaction is increased) and sufficiently 
weak inter-component interaction (i.e., 512 < 1). Finally, we examined the dynamical evolution of the instability 
of such interlaced structures, which typically resulted in the destruction of the waveforms, in favor of simpler, more 
stable dynamical patterns. 

Nevertheless, there is still a number of important open questions for future consideration. For instance, it would be 
particularly interesting to examine whether it would be possible for the inter-component coupling to actually stabilize 
structures that are dynamically unstable in the single-component setting. Also, it would be useful to possess a 
systematic classification of the solutions (interlaced and non-interlaced ones) available in the multi-component system 
setting, similarly to the one-component classifications of [s^ [35|. Such efforts are currently underway and will be 
reported in future publications. 
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